Model of Multiphoton Transitions in a Current-Biased Josephson Junction 
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We present a simple model to describe multiphoton transitions between the quasi bound states 
of a current-driven Josephson junction. The transitions are induced by applying an ac voltage with 
controllable frequency and amplitude across the junction. The voltage induces transitions across 
the junction when the frequency u) satisfies nhui = A_Eio, where A_Eio is the splitting between the 
ground and first excited quasi-bound state of the junction. We calculate the matrix elements of 
the transitions as a function of the dc bias current I, and the frequency u> and amplitude Vac of 
the microwave voltage, for representative junction parameters. We also calculate the frequency- 
dependent absorption coefficient by solving the relevant Bloch equations when the ac voltage is 
sufficiently weak. In this regime, the absorption coefficient is a sum of Lorentzian lines centered at 
the n-photon absorption frequency, of strength proportional to the squared matrix elements. For 
fixed Vac, the transition rate for an n-photon process usually decreases with increasing n. We also 
find a characteristic even-odd effect: the absorption coefficient typically increases with / for n even 
but decreases for n odd. Our results agree qualitatively with recent experiments. 

PACS numbers: 



I. INTRODUCTION 

A current-driven Josephson junction can exhibit clear 
experimental manifestations of quantized energy levels [lj, 
S 0, IS, IE IE 0]- For example, if the junction is biased so 
that the first minimum of the washboard potential holds 
only a few quasi-bound states, a suitable ac current can 
induce transitions between the lowest levels. Population 
of the excited levels leads to tunneling through the bar- 
rier and hence a voltage pulse which can readily be de- 
tected experimentally. In part because of this behavior, 
such junctions are promising candidates for use as qubits 
("phase qubits") in quantum computation[8|. 

Fig. [T] shows a schematic of the level structure for 
such a junction. Since the potential near a minimum 
is anharmonic, the level spacings are unequal, as indi- 
cated in the sketch. If A_Eio = E\ — Eq is the spac- 
ing between the ground and first excited state, and 
A£^2i = E2 — Ei is that between the first and second 
excited level, then AE21 < AE10. An ac perturbation of 
frequency lj = AEio/h will induce transitions between 
the lowest two states, but because of the anharmonic- 
ity, will not produce further excitations to the next level. 
Rabi oscillations between the lowest two levels have been 
demonstrated experimentally - that is, the population of 
the excited level oscillates with a frequency related to the 
amplitude of the ac pulse, if that perturbation is tuned 
to be in resonance with A_EioQ. 

Recently, Wallraff et al9] have demonstrated multi- 
photon transitions between the ground and first excited 
states of a current-biased Josephson junction. In this 
experiment, transitions were induced between the two 
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levels by subjecting the system to an ac voltage at a 
frequency AEio/(nh), with n > 1. Similar transitions 
have also been observed in flux qubits by Saito et aljlQl]. 
and in general in qubits based on anharmonic potentials 
containing more than two quasi-bound states [11(. These 
transitions correspond to absorption of n microwave pho- 
tons. Such transitions can be used to generate photons at 
multiples of the incident frequency. If the states involved 
have a sufficiently long coherence time, the perturbation 
could also lead to Rabi oscillations arising from multipho- 
ton transitions. Indeed, such Rabi oscillations have been 
reported experimentally [l~2| . Thus, they could be a sig- 
nificant advance in using these phase qubits in quantum 
computation. 

A number of workers have discussed models for Joseph- 
son junctions driven by strong ac perturbations. For 
example, Ashhab et a/. [l3[ have considered a two-level 
system subject to a strong perturbation which is har- 
monic in time, and demonstrate the occurrence of res- 
onances at frequencies 00 = Q/n, where u> is the drive 
frequency, Ml is the unperturbed level splitting, and n is 
a positive integer. Saito et aL[12l| have observed and an- 
alyzed Rabi oscillations corresponding to n-photon tran- 
sitions in a superconducting flux qubit. Inomata et q/[l4j 
have measured and analyzed macroscopic quantum tun- 
neling in intrinsic P^S^CaC^Os+a Josephson junctions 
arising from multiphoton transitions, using an analysis 
based on the classical Josephson junction equation of mo- 
tion. Koval et aL[15|| have analyzed the enhancement of 
Josephson phase diffusion by microwaves, starting from 
the classical equation for a resistively and capacitively 
shunted Josephson junction. They obtain a character- 
istic dependence of the n-photon absorption cocficicnt 
on the square of an n th order Bessel function, similar 
to that obtained by earlier workers [HI, 0~7j], and also re- 
sembling that found in the different analysis presented 
here. Fistul et q£.[l8| have analyzed the quantum es- 
cape of the phase in a strongly driven Josephson junc- 
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tion in the presence of both ac and dc bias curents, in 
good agreement with their own experiments. Goorden 
and Wilhelm[l^| analyzed nonlinear driving effects in a 
continuously driven solid-state qubit using the Bloch- 
Redfield equation, once again obtaining n-photon reso- 
nance effects. Gronbech- Jensen et aL[2Cj] have obtained 
multiphoton-like effects in Josephson junctions computa- 
tionally, by solving the classical equation of motion for a 
junction in the presence of both dc and ac currents. A 
number of other workers have considered effects of strong 
ac fields on Josephson junctions in various geometries, 
mostly in the context of their possible use as qubits for 
quantum computation [U [22|, [23|, [24[ • 

In this paper, we provide a simple model to describe 
this multiphoton absorption in Josephson phase qubits, 
using a somewhat different approach from those de- 
scribed above. In our model, we consider a Josephson 
junction in the presence of a dc driving current plus an 
ac voltage. The junction is assumed to have very little 
dissipation, so that the only terms in the Hamiltonian are 
the Josephson coupling and a capacitive energy. These 
two combine to produce the well-known pendulum-like 
Hamiltonian of the junction. The spacing between the 
levels of this system depends on the current bias. 

When the ac voltage is introduced, the gauge-invariant 
phase difference in the Hamiltonian includes a term aris- 
ing from that voltage. Use of a Bessel function expansion 
then shows that this extra term produces n-photon tran- 
sitions between the ground and first excited state of the 
current-biased junction. We calculate the n-photon tran- 
sition rate as a function of ac voltage amplitude and dc 
driving current, using the Fermi golden rule. Our model 
readily produces the transitions seen in experiments. 

The remainder of this paper is organized as follows. In 
the next section, we briefly describe our calculation of 
the energy eigenvalues and wave functions of the time- 
independent Hamiltonian. Following this, we describe 
the expansion which leads to all the n-photon terms in 
the Hamiltonian. Finally, we present numerical results 
for the n-photon transition rate, using the Fermi golden 
rule. A brief discussion follows in Section V. 



II. FORMALISM 
A. Current-Biased Josephson Junction 

1. Hamiltonian 

In the absence of resistive shunting, a capacitively 
shunted Josephson junction driven by a current / can 
be described by the Hamiltonian |25j H = (2en) 2 / (2C) — 
coscf> — 7^-(t>, where I c and / are the junction critical 
current and the applied current, (f> is the phase differ- 
ence across the junction, C is the junction capacitance, 
n is the Cooper pair number difference across the junc- 
tion, and 2e is the charge of a Cooper pair, n and <p are 
canonically conjugate operators, satisfying the commu- 



tation relation [n, <f>] — —i. This relation can be satisfied 
if we use the representation n = —id /dip for the number 
operator. In this representation, H becomes 



d 2 hi 

H= - Ec w^ C0S( 



hi 
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-E c 



J^+£/W>), (i) 



where Ec = 2e 2 /C is the charging energy of the junc- 
tion. Eq. ([I]) describes a phase "particle" moving in a 
"washboard potential," U(4>) = — 4-^ cos0 — |^</>, whose 
slope is controlled by /. In the experiments of Ref. 
Ec -C Ej = hl c /(2e). Although this regime corre- 
sponds to a relatively large junction, it must still be 
treated quantum-mechanically when the dissipation is 
very small. For their parameters, the observed quasi- 
bound state transitions occur when I / I c ~ 0.99. In our 
calculations, for calculational convenience, we have used 
a larger ratio of Ec/Ej. In order to model the rele- 
vant experiments, we solve both the time-independent 
and the time-dependent Schrodinger equations for the 
above Hamiltonian, as we now describe. 



2. Solution of the Time- Independent Schrodinger Equation 

For I < I c , U(4>) has local potential minima at (j> = 
<f>k = sin -1 (I/I c ) + 2nk, k = 0, ±1, ±2,.... For a given 
/, the potential near the first of these minima varies as 
U{4>) ~ U m in + K(I)((t> - o ) 2 /2, where U min = U{4> ) 
and K(I) is a current-dependent second derivative. If the 
potential were strictly harmonic near the minimum, the 
solutions of the time-independent Schrodinger equation 
would be harmonic oscillator eigenstates, with energies 
E n (T) = (n + ~)flu) c (I). These frequencies are readily 
shown to satisfy 



W C (J) = W C (0) 



I 2 



1/4 



(2) 



where w c (0) = [2e/ c /(SC)] 1//2 is the harmonic oscillator 
frequency at / = 0, usually called the junction plasma 
frequency. Because the potential is actually anharmonic, 
there are two corrections to the above expression for the 
energies of the eigenstates: the junction levels are not 
equally spaced as in a harmonic potential, and they are 
only quasi-bound states, because they can tunnel out of 
the well. 

It is convenient to work with a scaled Hamiltonian 



H 



H 
~E~j 



COS ( 



I 
X 



E c d 2 
Ej dej) 2 ' 



(3) 



The parameters appropriate for the experiments of Ref. 
8] are Ej - 6.8 x 1CT 14 ergs and E c ~ 8.3 x KT 20 ergs. 
In the present calculations, we generally use somewhat 
different parameters, as discussed below. 

To solve the time-independent Schrodinger equation 



(4) 
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for a given I/I c , we first find the local minimum c/)q — 
sin _1 (//J c ) and corresponding potential U m in = U(4>o) 
of the normalized washboard potential U(cj>) = coscj) — 
j-(f). We also find the closest local maximum 4> max and 

the corresponding relative potential maximum, U max = 
U{4> m ax)- For I <~ I c , it is readily shown that the barrier 
height AJ7 = U max - tj min oc [1 - I/I c f /2 ■ The number 
of quasi-bound states Nb(I) in the well near (f> — (f>o is 
approximately N b (I) ~ AU{I)/[hiu c (I)}. The lifetime of 
a given quasi-bound state in the well is also determined 
by this ratio. 

To find the quasi-bound states, we expand the solu- 
tions of eq. J4} in the set of harmonic oscillator states 
corresponding to the local minimum of U (c6) for the given 
1 /Ic [20| ■ These take the form 



M<t>) = 2-"/vr 1/2 



1 



1/4 



x exp 



(A0) 2 7>\ 



(5) 



2(A0) 2 J'"\ A0 
where H n (£) is a Hermite polynomial, and 



/ j2 \ -1/4 

(A</>) 2 = [2e 2 /hC][hC/(2eI c )} 1 / 2 f 1 - ^ 



(6) 



In this form, the orthonormal states are harmonic oscil- 
lator eigenstates of the Hamiltonian 



Hharm (-0 



Mr 



u min + 2e M 



j2\ V2 



(7) 



where we have used the relation {2]). 

The solution ^i((f>) of the time-independent 
Schrodinger equation is then expressed as 
^i{4>) = 5Z n o, n ij) n {(j)) . In matrix form, the Schrodinger 
equation becomes 



(8) 



where H mn — ip m ((f))Hip n ((f))d(j). In our calcula- 
tions, we include as many as the first 23 harmonic oscil- 
lator states in our matrix. As a check of our procedure, 
we have calculated the lowest three energy eigenvalues 
Eq, Ei, and E2, and corresponding wave functions, for 
I / I c = 0.99, approximately the dc current studied in Ref. 
Q, using their quoted values of Eq and Ej, and from 
these compute AE%o, AE21 and the ratio AE21/ AE\q. 
We find Al^io = 3.86 x 10~ 17 erg, quite close to the 
value of 4.51 x 10~ 17 erg measured in Ref. ||. The dif- 
ference may arise because the actual values of I/I c and 
C in the junction studied experimentally differ slightly 
from the quoted experimental values. We also calcu- 
late AE 2 i/AE w — 0.928, quite close to the value of 0.9 
quoted for this quantity in Ref. Q ■ 



Once we have these eigenstates, we determine how 
many of these are "quasi-bound" by finding out how 
many satisfy E n < U max . In general, we have consid- 
ered only values of I/I c such that there are at least three 
quasi-bound states. For the parameters of Ref. [8j, this 
condition is satisfied up to I/I c ~ 0.99. 



B. Current-Biased Josephson Junction with AC 
Voltage 

In the presence of an ac voltage V ac sinujt, t he p roper 
gauge-invariant phase difference takes the form[25]: 4> — > 
cos ujt. The correspondingly modified Hamilto- 



2eV a , 

nian is H 

V ac I 



-E, 



9 2 



hi, 

2e^ 



[4>- 



2c.V n , 



coswi 



; 90 2 2e^ 2e 1 huj 

■ cos tot. This formulation makes clear that the ef- 
fects of the ac voltage are characterized by the frequency 
u> and by the dimcnsionless variable z — 2eV ac /(huj). In 
the dimensionless form of Section IIA, the Hamiltonian 
takes the form 



H 



Ec_cP_ 
Ei d6 2 



I 

T 



> + — §^ cos(a)r) 



cos(d)r), (9) 



where we have introduced the dimensionless frequency 
q _ 2ew dimensionless time r = 4f^, and dimensionless 

l c ' 2e ' 

ac voltage amplitude V ac = j^j-, where Rq — fr/(4e 2 ). 

We are interested in terms in this Hamiltonian which 
will induce transitions between the quasi-bound states 
discussed in Section IIA. Since the last term in eq. 
(|9"|) is a c-number, it will not induce such transi- 
tions. The only relevant term is the Josephson en- 



ergy cos 



+ Sr cos(u)t) 



To extract the multipho- 

ton transitions, we may express this energy in terms of 
Bessel functions using the expansions (2?| cos(zcos(#)) = 
J (z) + 2Y l %Li(-) kj 2k(z)cos(2k6) and sin(zcos(6»)) = 
2 Er=o(-) fej 2fc+iWcos(2fc+l)0. Writing z = V ac /Q 
and 9 — ujt, and using a trigonometric identity, we can 
easily show that 



cos(4> + -z- cos(q)t)) = Jo(—^- I cos < 



oo 



+2 \^ J n I — — ) cos(0 + nn/2) cos(n(Dr). 

n—l 



(10) 



Thus, in the presence of an ac voltage, the potential en- 
ergy term in the Schrodinger equation is modified in two 
ways: (i) the strength of the dc part is reduced by a fac- 
tor of Jo(V ac /ui), and (ii) there are an infinite series of 
additional ac terms at frequencies uj and all its harmon- 
ics. These harmonic terms will induce the multiphoton 
transitions. Given these approximations, the full Hamil- 
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tonian takes the form H = Hq + Hi , where 



d 2 hi hl c 

~ Ec W~Te (t) -^ jQ 



2eK, 



COS0 (11) 



and 



Hi 



2eV ai 



cos(0 + nir/2) cos(nwt). 



The condition for the occurrence of an n-photon tran- 
sition between the ground and first excited state of the 
current-driven junction is 



AE 10 (I, V aci w) = nhw. 



(13) 



Here AEio(I,V ac ,u) is the energy splitting between the 
ground and first excited state when the ac voltage am- 
plitude and frequency are V ac and u>. We solve eq. (fT3|) 
by a self-consistent procedure. For a given V ac and u>, 
the junction levels are obtained as solutions of the time- 
independent Schrodinger equation H(V ac , w)^ = E^f, 
where H(V ac ,uj) is given by eq. (2) but with Ej re- 
placed by EjJ [2eV ac /(hui)] To find the strength of the 
ri-photon transitions, we make an initial guess for w, 
solve the dc Schrodinger equation using the potential 
— Ej[Jo(2eVac/(hw)] cos (f>, and iterate until eq. (fT3|) is 
satisfied. As in the previous section, we expand the quasi- 
bound eigenstates as a linear combination of harmonic 
oscillator wave functions. Thus, we make an initial guess 
for u>, calculate the eigenvalues, and then change u> such 
that nhw = AE 10 (I,V ac ,w). The eigenvalues are recal- 
culated using this new u>, and the procedure is repeated 
until the energy levels remain unchanged to within an 
absolute error of 1.0 x 1CT 4 Ej. Although eq. (fTU)) is 
valid, in principle, for any value of z, this procedure 
is not reasonable unless the strength of the potential 
[M c / (2e)]Jo[2eV ac / (hu>)} > 0. In practice, we find that, 
for our choice of parameters (see below) , our procedure 
leads to two or more bound states in the well only if 

2eV a , 



fru) 



< 0.58. 



Given the energy levels, the n-photon transition rate 
can be calculated, in principle, from the Fermi Golden 
Rule[H||: 

2ir 

T n {I,V ac ,w) = — \(0\V n \l)\ 2 5{E 8 -nhw) 
n 

= T n (I,V ac ,w)S(AE w -nhw), (14) 

where the perturbing potential for n-photon transitions 
is V n = 2J„(^y^)cos(</> + mr/2). The states (0| and 
|1) are the ground and first excited states of the current- 
biased junction in the presence of microwave radiation. 

Eq. (fT4|) is suitable for describing transitions between 
bound states with sharply defined energies. In reality, 
as already noted, the states of interest are only quasi- 
bound, since, for both the initial and the final state, the 
phase "particle" can tunnel out through the barrier of 
the washboard potential. One way to take account of 



this tunneling is to describe both the initial and the fi- 
nal states by complex energies E + ij and E\ + vyi. 
The quantities 70/ft and Ji/H then represent the tunnel- 
ing rates out of the states [0) and |1). If this description 
is used, the n-photon transition rate is described by the 
second form of the Fermi Golden Rule[28(, which repre- 
sents the transition rate from an initial state of energy 
Eq into a continuum of final states: 

T n (I,V ac ,w) = ^-\{0\V n \l)\ 2 p{AE 10 - nhw). (15) 

Here p(E) is the density of final states. If the final state 
is approximated by a complex energy Ei + 171 , then the 
quantity which replaces the delta function in eq. ()14l) is 



p(AEiQ — nhw) 



7i 



1 



ir (AE W - nhw) 2 + j 2 



(16) 



where the normalization is chosen so that p(E)dE = 
1. The corresponding rate of energy absorption is then 



a = y^ nhwT n (I, V ac ,w)(nHw), 



(17) 



where the factor of nhw denotes the fact that the energy 
absorbed in an n-photon transition is nhw. 

In order to use this formulation, we need to determine 
the width 71 . For the regime studied in typical experi- 
ments, it should be sufficient to treat 71 within the WKB 
approximation. This approximation gives 



71 ~ hwQ exp {—2S/K) , 



(18) 



where wq is a suitable attempt frequency for escap- 
ing the well, and S is the WKB action. In this case, 



S(E) = f 



4>f 



where — —ih(d/d<p) is the momen- 



tum canonically conjugate to (f>, and fa and fat are the 
values of 4> at the left and right hand edges of the tun- 
neling barrier. From the Hamiltonian (1) of the current- 
driven Josephson junction, we replace p$ in the WKB 
approximation by h^J\E — U{fa)]/Ec- Thus, we finally 
obtain (setting E = Ei for the final state) 



S(E) 



+' E - U{fa 
, V E c ' 



(19) 



where fa and far are determined by U(fa) = U(far) = Ei. 
We give some numerical estimates of 71 , and of the corre- 
sponding multiphoton absorption, at the end of Section 
IV. 

Besides the escape rate from the metastable poten- 
tial well, another source of line broadening is dissipation 
due to the finite shunt resistance R of the underdamped 
Josephson junction. We expect that this dissipation will 
give rise to another contribution to 71, thereby further 
increasing the linewidths of the multiphoton absorption 
lines. The magnitude of this contribution obviously de- 
pends on the magnitude of R. 
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III. INCLUSION OF DISSIPATION USING 
BLOCH EQUATIONS 

The form of the result (fT7|) for multiphoton absorption 
in the presence of dissipative processes can also be ob- 
tained using the analog of the Bloch equations familiar in 
NMR. In this section, we describe this method. Related 
approaches have been discussed, e. g., by Shevchenko et 
«l[29j] . Our approach goes beyond their work because we 
directly compute the necessary matrix elements and thus 
obtain explicit expressions for the absorption coefficient. 

We again consider the Hamiltonian H = Hq + Hi of 
cqs. (fTTj) and (fT2|) . and consider transitions only between 
the lowest two states of the Hamiltonian (fTTj) . Since we 
are considering only two states, we may write the Hamil- 
tonian (fTTj) in operator form as Hq = — A ^ 10 <J Z + const, 
where AEi (I, V ac ,uj) denotes the splitting between the 
lowest two energy levels in the presence of an ac voltage 
of amplitude V ac and frequency u>. The last term is a 
constant which can be chosen to vanish by proper selec- 
tion of the energy zero, we choose our zero of energy so 
that this constant vanishes. Ui (i = x, y, z) are the three 
standard Pauli matrices. If we retain only that part of 
H\ which produces transitions between |0) and |1), then 
the time-dependent perturbation (|12[) can be written 



Hi — ^ 2 J n 



2eV ac 



| cos(4> + nir/2)\l) cos(ncut)a x 



In writing the Hamiltonian in this form, we are not only 
considering just two levels but also are using the fact that 
the diagonal matrix elements of the form (0\Hi(t)\0) and 
(l|#iO)|l) vanish. 

In the absence of relaxation processes, the Heisen- 
berg equations of motion can be used to obtain equa- 
tions of motion for expectation values of any opera- 
tors. We write the time-dependent wave function as 
|^(t)) = a (t)\0) + b(t)\l), where \a\ 2 + \b\ 2 = 1. We 
also write the Hamiltonian as H = Hq + Hi = — B • er/2, 
where a is the ordered triple of Pauli matrices, and B is 
an effective magnetic field, which is the sum of a time- 
independent B z z and a time-dependent part B x (t)yi, with 



and By = 0. 



B 7 = A£i 



2eV n , 



o, (21) 

cos(0+7i7r/2)|l) cos(nwi), 
(22) 



The Heisenberg equations of motion for the ex- 
pectation value of any operator O are ih-^(0) = 
ih-^(ip\0\tp) = {ij)\[0,H]\ijj). For the present Hamilto- 
nian H = Hq + Hi , including only the two lowest energy 
levels, they reduce to the well-known form 



4-M = Im x B, 

dt h 



(23) 



where M = (a), the brackets denoting a quantum- 
mechanical expectation value. 

The form of eq. (|23p suggests that, to include dissipa- 
tive processes, we should generalize this equation so that 
it has the same form as the Bloch equations of magnetic 
resonance: 



dt 
d 
It 



Mi 



h 
l 



(M x B) 



(M z - Mq) 



(24) 



Here and r r are relaxation times analogous to Ti and 
T2 in magnetic resonance theory, = {M x , M y ), and 
Mqz is the equilibrium value of (a) to which M reverts 
when the time-dependent perturbation is turned off [see, 
e. g., Ref. 

Given the solution for M(i), we can calculate the ab- 
sorbed power by adapting a standard approach used in 
NMR calculations [3(j, namely, we solve for M x (t) in the 
limit when B x {t) is weak. In this regime, M x (t) can be 
obtained in closed form by superimposing the solutions 
arising from each separate frequency f2 of B x (t) in eq. 
([22)1. For small B x , M z ~ Mq through first order in B x 
Writing the three Bloch equations in component form, 
and also assuming that both M x and M y vary with time 
as exp(— iQt), we obtain 



itiM x 



-iClM,, 



1 M x 

— MyB Z X ~ 

1 



{M B X - M X B Z 



The resulting solution for M4 
found to be 



M+(t) 



(B x /h)M e 



(B z /h)-n-%/T4 



iM y is readily 
(26) 



For the present problem, the transverse field is a sum of 



terms of the form B nx cos(nwi) = (B nx /2)(e 
which induces a transverse magnetization 

M nt+ {t) = {B nx /2h)M Q x 

/ ^ — incut 

X 



')• 



B„ z /h - nuj - i/r^ B z /h + nuj - ijr^ 



(27) 



Hence M nx (t) = ReM„ !+ (i) satisfies 

2HM M n 



cos(nu;t)(B nx /h — ntu) + sm{nujt){l / r^) 

(B z /h-nLu) 2 + l/T 2 
cos(nujt)(B nx /h + nuj) — sin(nwi)(l/r0) 
+ (B z /h + nu;)i + l/T 2 ' 

Finally, the total absorption can be written 



a = (-(Hi(t))) t 



(28) 



(29) 
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where the outer brackets denote a time average. For our 
time-dependent perturbation, this expression is readily 
shown to be equivalent to 



^Mo , Bnx 

n=l 



(B z /h-nu;) 2 + l/T 2 (Bz/h+nu) 2 + 1/rJ 



(30) 



In terms of the parameters of the two-level system, we 
have 



where 



Mo 



2e 2 



7'Jn 



2eV ai 
tiu> 



nuj 
4r 



(31) 



|cos((^+n7r/2)|l)| 2 £( nuj), 



where C{nu>) is a Lorentzian given by 

1 

C{nio) = 
+ 



(32) 



(AE 10 /h - ncj) 2 + 1/r 2 , 
1 

(AE w /h + nLu) 2 + l/ T 2 ' 



(33) 



Thus, in this approximation, the total absorption is just 
a sum of Lorentzians. The quantity Mq, undetermined 
in this calculation, is given by 2p — 1, where p is the 
probability that the junction is to be found in the lower 
of the two states, |0). The integrated strength of the n- 
photon absorption line is 



f°° . . 2vr 2 

/ a„ w = Mq-t-t&B B nz = F T , 

J-oo n, 



(34) 



where T„ is the coefficient of eq. ((HI). Eqs. dSU, (021), and 
(l33|) are basically equivalent to the result (flTf obtained 
in the previous section. The quantity h/r^ is analogous 
to the quantity in eq. (JTSJ) . 

In principle, a more accurate solution could be ob- 
tained by solving the Bloch equations directly for the full 
time-dependent Hamiltonian arising from the microwave 
voltage. This would lead to a more complicated absorp- 
tion lineshape, which would probably depend on r r as 
well as T^. 



IV. NUMERICAL RESULTS 

We now describe our calculated n-photon transition 
rates T n (I,V ac ,u>), in the absence of dissipation, for n 
ranging from 1 to 5. In each calculation, we start by 
choosing /, V ac and n. The frequency and energy-level 
splitting are then determined self-consistently, as de- 
scribed above, and the transition rate T n is calculated 



[see eq. (|14p]. In all calculations, we choose Ec = 
8.53 x 1(T 18 erg, Ej = 6.91 x 1CT 14 erg, and hence 
Ec/Ej = 1.235 x 10~ 4 . For this choice, there are 
three quasi-bound states at I/I c ~ 0.85, and the ratio 
AE 2 i/AE w ~ 0.928. 

In Fig. [21 we show Y-y as a function of z = 
2eV ac / (hco) for a narrow range of driving currents, 0.85 < 
I/I c < 0.856. The frequency is chosen so that hoj — 
AEiQ(I,V ac ,tvjj). For the chosen range of /, Ty is quite 
insensitive to /. In Fig. [31 we plot AEia(I, V ac , cv) as a 
function of V ac at a frequency satisfying hw = AE±q. 
The splitting decreases roughly quadratically with in- 
creasing V ac - This behavior is expected from the Bessel- 
function dependence of the potential strength, since Jo (z) 
decreases approximately quadratically with z, resulting 
in a shallower potential well with increasing z for fixed 
I. 

In Figs. 2] and O we show corresponding results for 
two-photon transitions (n = 2), all other parameters be- 
ing kept the same as in Figs. [2] and [3) In these plots, 
we have considered smaller values of V ac than for n = 1. 
We do so because of the condition nfru — AEi (I, V ac , u) 
which must be satisfied for rt-photon absorption. Two- 
photon absorption thus occurs at roughly half the fre- 
quency of one-photon absorption, as expected. AE\$ is 
more sensitive to V ac at such frequencies than it would be 
for one-photon absorption. We must also keep V ac well 
below the limit V ac = 0.58nhuj/ (2e), above which there 
are no quasi-bound states. 

The corresponding results for three-photon transitions 
are shown in Figs. [S] and [7] In general T3 < Ti, as ex- 
pected, since the corresponding matrix element involves 
Js{z) instead of J\(z). Furthermore, of course, for a given 
I and Vac, the frequency required to produce the n = 3 
transition is smaller than that for n = 1. We have also 
considered smaller values of V ac than for n = 1, because 
the sensitivity of AE\q to V ac is greater for n = 3 than for 
n = 1. In Figs. 181 -[TT1 we show the corresponding results 
for n = 4 and n = 5 transitions. 

All the above results make clear a noticeable differ- 
ence between even and odd transitions: for n even, T n 
increases with increasing / at fixed z, whereas for n odd, 
r„ decreases with increasing /. (We have not considered 
a very large range of I in calculating these rates.) This 
difference is easily understood from the matrix elements 
of the time-dependent perturbation which enter equation 
(TT4"]) for the transition rate. Specifically, the argument of 
the cosine term has a phase shift of nrr/2, so that the per- 
turbing potential varies as sin <f> for odd n. Because the 
perturbation is even in (f> for even n, T n vanishes for even 
n at I = 0, and hence should increase within increasing 
/, as we observe. By contrast, the perturbation is odd 
in <f> for odd n, so that T n for odd n is finite at I = 
but decreases with increasing /. Thus, we interpret this 
even-odd effect as arising from the different parities of 
the time-dependent perturbation for n even and n odd. 

Figs. 2-10 also show that, for fixed z = 2eV ac / (hu>) 
and fixed I/I c , r„ generally decreases with increasing 
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n. Comparison of the behavior for different n is some- 
what difficult, because the relevant value of the parame- 
ter z = 2eV ac / {hw) is different for each n. This behavior 
is consistent with that observed in the experiments of 
Ref. [8J. These authors find that, when the first excited 
state is populated by n-photon transitions, the appar- 
ent lifetime increases with increasing n. This increased 
lifetime is simply the consequence of a smaller transition 
rate for the higher-n transitions. 

Finally, we estimate how these n-photon transition 
rates are affected when we include reasonable estimates of 
level broadening due to tunneling of the phase "particle" 
through the barrier. We use the formalism of eqs. (fT5|) - 
(HHJ). Using the values of Ec and Ej given above (which 
lead to Eq/Ej = 1.235 x 10~ 4 ), we obtain estimates for 
exp[-2S(£a)/ft] given in Table 1 for I/I c = 0.85, 0.90, 
and 0.95. The corresponding width 71 of the Lorentzian 
line shape is given, as a fraction of hio c (I), in the last 
column. In calculating S(Ei), we have estimated E\, 
the energy of the second-lowest quasibound state, by ex- 
panding the potential in the Hamiltonian (1) through 
third order in cf> — <j>o, and have included the cubic term 
as a perturbation to the resulting harmonic Hamiltonian, 
through second order in perturbation theory. The result- 
ing values of E\ are also shown in the Table. 

In order to illustrate a typical multiphoton absorption 
spectrum in the presence of this dissipation, we have car- 
ried out a simplified model calculation of the total ab- 
sorption spectrum, f tat (I, V ac , uj) = J2nZT ^n(I, Vac, w), 
where T n is given by eqs. (fTS"|) and (fl6|) . The results 
are shown in Figs. 12 and Fig. 13 for I/I c = 0.95, 
E c = 8.53 x 10" 18 erg, and Ej = 6.91 x 10~ 14 erg as 
earlier in this paper. We assumed I/I c = 0.95 and con- 
sidered two ac amplitudes given by 2eV ac /hwo(I) = 0.2 
and 0.4. We have included all multiphoton terms through 
n ma x — 7. For this calculation, we have calculated both 
the energies Eo and E\ and the matrix elements (0|V^|1) 
by including only the lowest order anharmonic correc- 
tions to Eo, Ei, |0) and |1), the perturbation being the 
terms in the potential of third order in <fi — 4>q. In the 
interests of simplifying this calculation, we have also ne- 
glected the dependence of the Josephson well depth on 
the ratio z — 2eV ac /{huj) [eq. (11)]. In practice, this 
means that this model calculation becomes increasingly 
inaccurate at low frequencies. For our choice of V ac , 
this dependence on V ac should probably be included for 
hw < 0.2A£io. With these assumptions, however, both 
the energies and the matrix elements can be computed 
in closed form in terms of Ej, Ec and I/I c . The results 
show that the strength of the lines falls off rapidly with 
increasing n for our chosen value of V ac . 

In these calculations, and with our estimate of 71, we 
have not included the effects of the shunt resistance R, 
which is always present in a realistic Josephson junction. 
This shunt resistance will further broaden the multipho- 
ton absorption spectrum, but is distinct from the tunnel- 
ing through the barrier treated above. For the junctions 
studied in Ref. 0, l/[w p (0)r] < 5 x 10 -4 , using their 



parameters of critical current density j — 1.1 kA/cm 2 , 
junction area ~ (5.5/j,) 2 , C ~ 1.6 pF, and shunt resis- 
tance R > 500S1. Here r = RC is the relaxation time 
associated with resistive dissipation within the junction. 
If we assume the same value of 1/[w p (0)t] for our model 
junctions, the broadening of the multiphoton lines at 
Ijl c = 0.95 would be about ten times smaller due to this 
source of dissipation than that due to tunneling through 
the barrier. For other junctions, of course, this might be 
the dominant source of broadening. 



V. DISCUSSION 

We have described a simple model for multiphoton 
transitions in a current-biased Josephson junction, in the 
presence of microwave irradiation. The model allows 
calculation of the transition rates between the ground 
and excited states as a function of dc bias current, ac 
voltage amplitude, and frequency. In the absence of 
damping, the absorption occurs as a series of delta- 
function lines satisfying the energy-conservation require- 
ment nfojj = AEio(I, V ac , ui). The lines are broadened by 
dissipative processes. We calculate this broadening using 
the analog of the Bloch equations as applied to the lowest 
two levels of the junction. This directly gives the absorp- 
tion coefficient ol(uj), which is approximately a sum of 
Lorentzian absorption lines centered around frequencies 
corresponding to n-photon absorption. 

Experimentally, multiphoton transitions are detected 
via enhanced tunneling from the excited quasi-bound 
states through the barrier of the washboard potential into 
the continuum 8j. In the present work, we calculate the 
steady-state absorption rate a(uj), using the Bloch equa- 
tions. By conservation of energy, we expect that this rate 
should equal the rate of tunneling through the barrier 
from the state E±. Since our calculated a n {uS) generally 
decreases with increasing n, this would imply that the 
rate of tunneling through the barrier would also decrease 
with increasing n, as is reported in experiments Q. How- 
ever, we have not attempted to calculate the relevant 
phase relaxation time which determines the width of 
the Lorentzian peaks in this approximation. 

The calculation of a(u>) via the Bloch equations 
is presumably derivable from the master-equation 
approach [31| previously used to discuss the time- 
dependence of level occupation numbers in Josephson 
junctions. In the present work, we also explicitly cal- 
culate the matrix elements needed to calculate a n (u>), 
which allows us to explicitly estimate the multiphoton 
absorption rate as a function of /, V ac , and ui. 

Our calculations have several characteristic qualitative 
features. For example, in agreement with experiment Q, 
the n-photon absorption rate, for a given / and V ac , gen- 
erally decreases with increasing n. Also, if the damping 
is weak, the absorption will be small unless AE W ~ nHu. 
Finally, the integrated strength of the n-photon transi- 
tion exhibits an conspicuous odd-even effect: the transi- 
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tion rate from the ground to the first excited state gener- 
ally increases with increasing / for even n, but decreases 
for n odd. Indeed, at zero bias current, the transition rate 
vanishes for even n. This is just a consequence of the sym- 
metry of the eigenstates in the well: if the bias current is 
zero, the ground and first excited states have even or odd 
parity, thereby allowing only odd-photon transitions. 

In future work, one calculation of interest would be 
to solve the Bloch equations directly, i.e., without ex- 
panding the time-dependent potential in a Bessel func- 
tion series, and without making the small-amplitude ap- 
proximation which allows the absorption line to be de- 
composed into a sum of Lorentzians. Such a calculation 



should be straightforward to do numerically. Finally, of 
course, it would also be valuable to explicitly compute the 
additional contribution to the broadening arising from 
the shunt resistance in the junction. 
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FIG. 1: Schematic diagram of the level structure for a Joseph- 
son junction subjected to an applied dc current. 
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i 


01 


4>2 




ji/[tkj (l)} 


Ei/Ej 


0.85 


1.2632 


2.6636 


6.584 


3.04 x 10" 7 


-1.3762 


0.90 


1.3557 


2.4454 


4.740 


1.22 x 10~ 5 


-1.4336 


0.95 


1.5236 


2.1451 


2.349 


2.35 x 10" 3 


-1.4946 



TABLE 1: WKB action SiE^/h for three different values of 
the applied current i — I/I c , as indicated. The second and 
third columns of the table denote the values of the phase 4> 
at the left and right hand edges of the barrier at energy E\ . 
E\ is the energy of the second quasi-bound state in the well, 
as approximated by the method described in the text. The 
fourth column is the approximate ratio of the linewidth to 
the current-dependent small-oscillation frequency uio(I). The 
last column gives the energy of the first excited state Ei(I), 
in units of E.j, including anharmonic corrections to lowest 
non- vanishing order in <f> — <f>o . 
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FIG. 2: (Color online.) Single photon transition rate Fi (in 



units of Ej/h) as a function of z 



for currents ranging 



from I/I c = 0.85 to 0.856. The red (dark gray full) line, 
green (light gray dashed) line, blue (gray dotted) line, and 
violet (gray dash-dotted) line correspond to z = 0.85, 0.852, 
0.854, and 0.856 respectively. The same color and grayscale 
codes are used in Figs. 3-11. 
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FIG. 3: (Color online.) Energy level splitting AiJio (in units 
of Ej )between the ground and first excited states as a func- 
tion of Vac for currents I/I c ranging from 0.85 to 0.856 for 
the self-consistently determined frequency ui = AEio/h. Pa- 
rameters are the same as in Fig. [2] 
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FIG. 4: (Color online.) Two photon transition rate T2 plot- 
ted as a function of z = 2e ^° for currents ranging from 
I / I c — 0.85 to 0.856 and the self-consistently determined fre- 
quency uj = AEio / (2h) . 
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FIG. 5: (Color online.) Energy level splitting AEio as a 
function of V ac for currents ranging I /I c — 0.85 to 0.856 and 
a frequency such that uj — AE\o/(2h). 





FIG. 7: (Color online.) Same as Fig. [5]but for three-photon 
transitions. 
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FIG. 8: (Color online.) Same as Fig. [4] but for four-photon 
transitions 
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FIG. 9: (Color online.) Same as Fig. [5] but for four-photon 
transitions. 
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FIG. 10: (Color online.) Same as Fig. [4] but for five-photon 
transitions. 
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FIG. 11: (Color online.) Same as Fig. [5] but for five-photon 
transitions. 
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FIG. 12: Plot of the absorption spectrum f tot{I, V ac ,u) = 
YZ=T Tn(I, V ac ,u) (in units of AE 10 /h) as a function of w, 
for I/I c — 0.95 and ac amplitude given by 2eV a c/{AEio) = 
0.2, as calculated using the parameters and method described 
in the text. 




FIG. 13: Same as Fig. 12, but with ac amplitude given by 
2eV ac /(AE w ) = 0.4. 



